function [sx, sy, theta, H, CS] = RANSACaffine(X1, X2)
% set RANSAC options
options.epsilon = 1e-6;
options.P_inlier = 1-1e-4;
options.sigma = 0.1;
options.est_fun = @estimate_affine;
options.man_fun = @error_affine;
options.mode = 'MSAC';
options.Ps = [];
options.notify_iters = [];
options.min_iters = 1000;
options.max_iters = 10000;
options.fix_seed = false;
options.reestimate = true;
options.stabilize = false;

X = [X1; X2];
[results, options] = RANSAC(X, options);% run RANSAC
if numel(results.Theta) == 0
    sx = NaN;
    sy = NaN;
    theta = NaN;
    H = [];
    CS = [];
    return;
end
R = [results.Theta(1) results.Theta(3); ...
    results.Theta(2) results.Theta(4)]; %rotate and scale
T = results.Theta(5:6); % shift

% X2 should be R*X1+T
H = [R T;0 0 1];
sx = norm([R(1,1),R(2,1)],2);
sy = norm([R(1,2),R(2,2)],2);
theta = asin(R(1,2)/(sy + eps));
CS = results.CS;
end
